library(tidyverse)
library(terra)
library(sf)
library(rnaturalearth)
library(here)
library(iucnredlist) ### remotes::install_github("IUCN-UK/iucnredlist")
library(janitor)
library(cowplot)Species change at pixel level
See how species change at the pixel level
Will create a map of Soresen’s Similiarity Index at the pixel level
# Functions
source(here("common_fxns.R"))# Read in Lawler data
spp_mapfiles <- list.files(here_anx('/data_lawler_2020/species_projections'),
full.names = TRUE, recursive = TRUE, pattern = '.tif')
spp_all_df <- data.frame(f = spp_mapfiles,
spp = basename(spp_mapfiles) %>%
str_remove('(_gf85|_in85|_mc85|_bias).+') %>%
str_replace_all('_', ' '),
tax = basename(dirname(spp_mapfiles)))# Read in AVONET dataset (reading in BirdLife tab, which should match with IUCN names)
avonet <- read_csv(here("data-laa", "AVONET-birdlife-tab-supplementary-dataset-1.csv")) %>%
clean_names()# Filter for just Lawler birds
birds_all_df <- spp_all_df %>%
filter(tax == "Birds")birds_all <- birds_all_df %>%
group_by(spp, tax) %>%
summarize(files = n())Gather maps and aggregate
Following the method in 2_combine_spp_scenarios.qmd: For all the Lawler birds species (301) gather all future scenario maps and aggregate, summing and dividing by 3 to make a “probability” map for future distributions, and crop to just the lower 48 US states (and a little Canada and Mexico). Then crop the historical maps to the same extent.
Aggregate future projections
For each species in the sample, aggregate the three future scenarios into a single map. Crop to the lower continental US states.
### set up extent for continental US
us_extent <- ext(c(xmin = -130, xmax = -65, ymin = 23, ymax = 50))
### Aggregate over the vector of species names
spp_vec <- birds_all_df$spp %>% unique()
future_rast <- lapply(spp_vec, FUN = agg_scenario, df = birds_all_df) %>%
setNames(spp_vec) %>%
rast()
x <- future_rast[[1]]
#plot(x, main = names(future_rast)[1])crop historical maps
hist_rast <- lapply(spp_vec, FUN = agg_scenario, df = birds_all_df,
scenario = 'historical') %>%
setNames(spp_vec) %>%
rast()
x <- hist_rast[[1]]
#plot(x, main = names(hist_rast)[1])Generate jurisdiction map
Let’s generate a simple jurisdiction map - US states, and Canada and Mexico.
states_sf <- ne_states(returnclass = 'sf', country = c('Mexico', 'Canada', 'United States of America'))
state_id_df <- states_sf %>%
st_drop_geometry() %>%
mutate(juris = ifelse(adm0_a3 == 'USA', name, admin)) %>%
dplyr::select(id = gn_id, juris)
### plot(states_sf %>% select(iso_a2))
states_rast <- rasterize(states_sf, hist_rast, field = 'gn_id')
#plot(states_rast)# Function to get pixels instead of jurisdictions
extract_jurisdictions_pixels <- function(spp, future_r, hist_r, juris_r) {
### spp <- spp_vec[1]
### future_r = future_rast; hist_r = hist_rast; juris_r = states_rast
s_f_r <- future_r[[names(future_r) == spp]]
s_h_r <- hist_r[[names(hist_r) == spp]]
s_df <- data.frame(values(s_h_r),
values(s_f_r),
values(juris_r)) %>%
setNames(c('historical', 'future', 'id')) %>%
mutate(pixel_id = row_number()) %>%
filter(!is.na(id)) %>%
filter(!is.na(historical) | !is.na(future)) #%>%
#group_by(id) %>%
#summarize(n_historical = sum(historical > .25),
#n_future = sum(future > .5))
return(s_df)
}Rough pass at classification?
For each spp let’s create a dataframe of jurisdictions and cell values…
state_spp_list <- lapply(spp_vec, FUN = extract_jurisdictions_pixels,
future_r = future_rast,
hist_r = hist_rast,
juris_r = states_rast) %>%
setNames(spp_vec) %>%
bind_rows(.id = 'spp') %>%
left_join(state_id_df, by = 'id') #%>%
### note, Canada and Mexico have multiple jurisdictions (provinces and states)
### so sum up all those instances
#group_by(spp, juris) %>%
#summarize(n_historical = sum(n_historical),
#n_future = sum(n_future),
#.groups = 'drop')Now that we have each species at the pixel level (pixel id), how many species are in each pixel for historical and future?
pixel_spp <- state_spp_list %>%
group_by(pixel_id, juris) %>%
summarize(n_spp_historical = sum(historical > .25),
n_spp_future = sum(future > .5)) %>%
mutate(n_spp_diff = n_spp_future - n_spp_historical)# Vast majority of pixels are losing total # of species (simply from numbers perspective)
ggplot(data = pixel_spp) +
geom_histogram(aes(x = n_spp_diff),
bins = 50) +
labs(title = "Bird species gain from historic to future by pixel",
x = "species change")Next step - sorensen similarity index (test a few pixels)
# Test out with a single pixel (141325)
pxl_141325 <- state_spp_list %>%
filter(pixel_id == 141325) %>%
# Assign future pixels > 0.5 to spp prescence
mutate(future = case_when(future > 0.5 ~ 1,
TRUE ~ 0))pxl_si <- pxl_141325 %>%
group_by(pixel_id) %>%
summarize(shared = sum(historical == 1 & future == 1),
historical_spp = sum(historical),
future_spp = sum(future),
ssi = (2 *shared)/(historical_spp + future_spp))# Try all pixels
pxl_ssi <- state_spp_list %>%
# Assign future pixels > 0.5 to spp prescence
mutate(future = case_when(future > 0.5 ~ 1,
TRUE ~ 0)) %>%
group_by(pixel_id, juris) %>%
summarize(shared = sum(historical == 1 & future == 1),
historical_spp = sum(historical),
future_spp = sum(future),
ssi = (2 *shared)/(historical_spp + future_spp))# Distribution of Soresen Similiarity Index for bird species
# Mean is 0.40 (moderate similarity between historic and future)
ggplot(data = pxl_ssi) +
geom_histogram(aes(x = ssi),
bins = 50) +
geom_vline(xintercept = mean(pxl_ssi$ssi)) +
annotate("text", label = "mean = 0.40", x = 0.30, y = 15000) +
#facet_wrap(~juris) +
labs(title = "Soresen Similarity Index btwn historic & future by pixel (n = 152,594)",
x = "SSI")# SSI stats by jurisdiction
ssi_juris <- pxl_ssi %>%
group_by(juris) %>%
summarize(mean_ssi = mean(ssi),
median_ssi = median(ssi),
max_ssi = max(ssi),
min_ssi = min(ssi)) %>%
ungroup()ssi_juris %>%
mutate(juris = fct_reorder(juris, mean_ssi, .desc = TRUE)) %>%
ggplot() +
geom_col(aes(x = juris, y = mean_ssi)) +
theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1)) +
labs(y = "mean pixel SSI",
x = "jurisdiction",
title = "Bird species mean pixel SSI by jurisdiction")Create sf of raster ids to map the SSI by pixel
# Create a blank raster with same extent as x = hist_rast[[1]]
id_rast <- rast(x)# Assign IDs and then convert to sf
values(id_rast ) <- 1:ncell(x)id_rast_sf <- id_rast %>%
stars::st_as_stars() %>%
st_as_sf() %>%
rename("pixel_id" = "Acrocephalus familiaris")# Merge in the SSI
id_rast_sf_ssi <- inner_join(id_rast_sf, pxl_ssi)ggplot() +
geom_sf(data = id_rast_sf_ssi %>%
filter(juris == "California"), aes(geometry = geometry, fill = ssi, color = ssi), size = 0) +
scale_fill_viridis_c(limits = c(0,1)) +
scale_color_viridis_c(limits = c(0,1)) +
labs(title = "Birds in CA Sorenen's Similiarity Index")ggplot() +
geom_sf(data = id_rast_sf_ssi, aes(geometry = geometry, color = ssi, fill = ssi)) +
scale_color_viridis_c() +
scale_fill_viridis_c() +
labs(title = "Birds Sorenen's Similiarity Index")In general, SSI values are as follows (https://www.statology.org/how-to-calculate-and-interpret-the-sorensen-dice-index/):
- 0.80-1.00: Very high similarity (communities share most species)
- 0.60-0.79: High similarity (substantial species overlap)
- 0.40-0.59: Moderate similarity (some shared species)
- 0.20-0.39: Low similarity (few shared species)
- 0.00-0.19: Very low similarity (minimal species overlap)
historic_spp_num <- ggplot() +
geom_sf(data = id_rast_sf_ssi, aes(geometry = geometry, fill = historical_spp, color = historical_spp)) +
scale_fill_viridis_c(option = "plasma", limits = c(1,180)) +
scale_color_viridis_c(option = "plasma", limits = c(1,180)) +
geom_sf(data = id_rast_sf_ssi %>% filter(historical_spp == 0), aes(geometry = geometry), fill = "gray40", color = "grey40") +
labs(title = "Historic number of bird species",
subtitle = "(0 species in gray)",
fill = "count",
color = "count")future_spp_num <- ggplot() +
geom_sf(data = id_rast_sf_ssi, aes(geometry = geometry, fill = future_spp, color = future_spp)) +
scale_fill_viridis_c(option = "plasma", limits = c(1,180)) +
scale_color_viridis_c(option = "plasma", limits = c(1,180)) +
geom_sf(data = id_rast_sf_ssi %>% filter(future_spp == 0), aes(geometry = geometry), fill = "gray40", color = "grey40") +
labs(title = "Future number of bird species",
subtitle = "(0 species in gray)",
fill = "count",
color = "count")plot_grid(historic_spp_num + theme(legend.position = 'bottom'),
future_spp_num + theme(legend.position = 'bottom'))ggplot() +
geom_sf(data = id_rast_sf_ssi %>%
filter(juris == "Canada"), aes(geometry = geometry, fill = ssi, color = ssi), size = 0) +
scale_fill_viridis_c(limits = c(0,1)) +
scale_color_viridis_c(limits = c(0,1)) +
labs(title = "Birds in Canada Sorenen's Similiarity Index")Add in Avonet data and investigate trait change at pixel level
Let’s start with the birds that easily merge with the Avonet dataset (270)
# See how well the Avonet trait data merges in
# 270 of 301 have matches - could work with these species for now since they should be able to match with IUCN
birds_avonet <- left_join(birds_all, (avonet %>% rename(spp = species1)))
lawler_avvonet_birds <- birds_avonet %>%
drop_na(sequence)
lawler_avvonet_birds_vec <- lawler_avvonet_birds$sppNow, will add in the bird species that did not merge as easily for a total of 299
# Read in crosswalk of lawler to avonet species (created with input from gemini)
lawler_avonet_xwlk <- read_csv(here("data-laa", "lawler_avonet_xwlk.csv")) %>%
clean_names() %>%
dplyr::select(lawler_spp:name_status) %>%
drop_na()# Match the lawler avonet xlk with trait data
# The two extinct species drop out
lawler_avonet_traits_missing <- inner_join(lawler_avonet_xwlk,
avonet %>% rename(avonet_spp = species1))# Bind the easily matched birds and the birds that were missing in the merge
lawler_avonet_all <- rbind(
lawler_avvonet_birds %>% dplyr::select(-tax, -files),
lawler_avonet_traits_missing %>% dplyr::select(-avonet_spp, -name_status)
)pxl_spp_list_trait <- inner_join(state_spp_list, lawler_avonet_all) %>%
filter(!is.na(id)) %>%
filter(!is.na(historical) | !is.na(future))Create map
us_spp_trait <- pxl_spp_list_trait %>%
dplyr::select(spp:juris, habitat:primary_lifestyle) %>%
filter(historical > 0.25 | future > 0.5) %>%
mutate(future = case_when(future > 0.5 ~ 1,
TRUE ~ 0)) %>%
mutate(habitat_density = case_when(habitat_density == 1 ~ "Dense",
habitat_density == 2 ~ "Semi-open",
habitat_density == 3 ~ "Open")) %>%
mutate(migration = case_when(migration == 1 ~ "Sedentary",
migration == 2 ~ "Partial",
migration == 3 ~ "Migratory"))# Combine the categorical columns into one group column to create unique groups
# Count up these trait groups by pixel
us_spp_trait_groups <- us_spp_trait %>%
unite(trait_group, migration, habitat_density, habitat, trophic_level, trophic_niche, primary_lifestyle) %>%
group_by(pixel_id, juris, trait_group) %>%
summarize(historical_trait_group_count = sum(historical),
future_trait_group_count = sum(future)) %>%
ungroup()# If we apply the sorensen index to traits, then we care about the presence / absence of a trait group and not the # of species
us_spp_trait_groups_ssi <- us_spp_trait_groups %>%
mutate(historical_trait_group_presence = case_when(historical_trait_group_count == 0 ~ 0,
historical_trait_group_count != 0 ~ 1),
future_trait_group_presence = case_when(future_trait_group_count == 0 ~ 0,
future_trait_group_count != 0 ~ 1)) %>%
group_by(pixel_id, juris) %>%
summarize(shared = sum(historical_trait_group_presence == 1 & future_trait_group_presence == 1),
historical_trait_groups= sum(historical_trait_group_presence),
future_trait_groups = sum(future_trait_group_presence),
trait_group_ssi = (2 *shared)/(historical_trait_groups + future_trait_groups)) %>%
ungroup()# Merge in pixel geography for map
id_rast_sf_us_ssi <- inner_join(id_rast_sf, us_spp_trait_groups_ssi)# Create map (299 bird species)
ggplot() +
geom_sf(data = id_rast_sf_us_ssi, aes(geometry = geometry, fill = trait_group_ssi, color = trait_group_ssi), size = 0) +
scale_fill_viridis_c(limits = c(0,1)) +
scale_color_viridis_c(limits = c(0,1)) +
labs(title = "Bird trait groups Sorenen's Similiarity Index")What’s the change from species SSI to trait SSI
# Filter the species dataframe down to on the avonet birds
pxl_ssi_avonet_spp <- state_spp_list %>%
# Assign future pixels > 0.5 to spp prescence
filter(spp %in% lawler_avvonet_birds_vec) %>%
mutate(future = case_when(future > 0.5 ~ 1,
TRUE ~ 0)) %>%
group_by(pixel_id, juris) %>%
summarize(shared = sum(historical == 1 & future == 1),
historical_spp = sum(historical),
future_spp = sum(future),
ssi = (2 *shared)/(historical_spp + future_spp))# Merge in the SSI
id_rast_sf_ssi_avonet_spp <- inner_join(id_rast_sf, pxl_ssi)Would expect trait SSI to generally to be greater than species SSI
pxl_ssi_dif <- left_join(us_spp_trait_groups_ssi %>%
dplyr::select(pixel_id, trait_group_ssi),
pxl_ssi %>%
dplyr::select(pixel_id, ssi)) %>%
mutate(ssi_diff = trait_group_ssi - ssi)# Merge in the pixel geography
id_rast_sf_us_ssi_diff <- inner_join(id_rast_sf, pxl_ssi_dif)ggplot() +
geom_sf(data = id_rast_sf_us_ssi_diff, aes(geometry = geometry, fill = ssi_diff, color = ssi_diff), size = 0) +
scale_color_distiller(palette = "BrBG", direction = 1, limits = c(-0.4, 0.4)) +
scale_fill_distiller(palette = "BrBG", direction = 1, limits = c(-0.4, 0.4)) +
# Show areas 0 species and trait overlap in pink
geom_sf(data = id_rast_sf_us_ssi_diff %>%
filter(ssi == 0 & trait_group_ssi == 0),
aes(geometry = geometry), fill = "pink", color = "pink") +
# Show areas with full species and trait overlap in navy
geom_sf(data = id_rast_sf_us_ssi_diff %>%
filter(ssi == 1 & trait_group_ssi == 1),
aes(geometry = geometry), fill = "navy", color = "navy") +
labs(title = "Trait SSI - Species SSI",
subtitle = "(pink is 0 trait and 0 species SSI; navy is 1 trait and 1 species SSI)")#Let’slook closer at the pixels where the difference is 0 Zeros may be areas where there are 0 species SSI and 0 trait SSI (no similarity btwn historic and future), 1 species SSI and 1 trait SSI (complete similarity btwn historic and future), or nonzero values.
pxl_ssi_diff_0 <- pxl_ssi_dif %>%
filter(ssi_diff == 0)
pxl_ssi_diff_0_summary <- pxl_ssi_diff_0 %>%
mutate(diff_type = case_when(trait_group_ssi == 0 & ssi == 0 ~ "zero",
trait_group_ssi == 1 & ssi == 1 ~ "one",
TRUE ~ "non-zero")) %>%
group_by(diff_type) %>%
summarize(count = n()) %>%
ungroup()
zero_trait_species_ssi <- pxl_ssi_diff_0 %>%
filter(trait_group_ssi == 0 & ssi == 0 )
zero_trait_species_ssi_pxls <- zero_trait_species_ssi$pixel_id# Look at the zero ssi pixels
zero_species_ssi <-state_spp_list %>%
# Assign future pixels > 0.5 to spp prescence
# filter(spp %in% lawler_avvonet_birds_vec) %>%
filter(pixel_id %in% zero_trait_species_ssi_pxls) %>%
mutate(future = case_when(future > 0.5 ~ 1,
TRUE ~ 0)) %>%
# Assigning presence / absence for these birds
filter(historical == 1 | future == 1)Plotting a few of the North Dakota birds
# Seems like it's one of the few birds showing up in future North Dakota (even though it seems unlikely to be there based on research)
plot(hist_rast$`Xenus cinereus`, main = "Xenus cinereus - historical")plot(future_rast$`Xenus cinereus`, main = "Xenus cinereus - future")plot(hist_rast$`Parus hudsonicus`, main = "Parus hudsonicus - historical")plot(future_rast$`Parus hudsonicus`, main = "Parus hudsonicus - future")plot(hist_rast$`Anser rossii`, main = "Anser rossii - historical")plot(future_rast$`Anser rossii`, main = "Anser rossii - future")plot(hist_rast$`Tympanuchus pallidicinctus`, main = "Tympanuchus pallidicinctus - historical")plot(future_rast$`Tympanuchus pallidicinctus`, main = "Tympanuchus pallidicinctus - future")